몬테카를로 방법

통계물리학을 비롯한 다양한 분야에서 몬테카를로(Monte Carlo) 방법을 이용한 시늉내기(simulation)가 널리 쓰이고 있다. 넓은 의미의 몬테카를로 방법은 매번 다른 결과가 얻어지도록 마구잡이 수를 이용해 결과를 얻는 것을 뜻한다. 간단한 예를 하나 들어보자. 종이 위에 반지름이 1인 원을 그려놓고, 원이 내접한 정사각형을 하나 생각하자. 임의로 –1에서 1사이의 값 [eqn: x]를 하나 택하고, 마찬가지로 –1에서 1사이의 또 다른 값 [eqn: y]를 택해서 점 ([eqn: x,y])를 얻자. 만약 [eqn: x^2 + y^2 < 1]를 만족하면, 이 점은 원의 내부에 놓여있다는 뜻이고, 이 부등식을 만족하지 못한다면 원의 외부에 점이 놓인다는 뜻이다. 위의 과정을 [eqn: N]번 반복했더니, 마구잡이로 택한 점이 원의 내부에 놓인 경우가 [eqn: N_1]번 일어났다면, 이 두 수의 비 [eqn: N_1 / N]는 다름 아닌 원의 면적과 정사각형의 면적의 비율이 된다. 즉, [eqn: N]이 늘어날수록 [eqn: N_1 / N]의 값은 [eqn: pi /4]로 접근하게 된다. 간단히 이 과정을 구현해 원주율 [eqn: pi]를 계산하는 코드를 적어 보았다.

import numpy as np, matplotlib.pyplot as plt

for N in 10**np.array([1,2,3,4,5,6,7]):

x = np.random.rand(N)*2 - 1.0

y = np.random.rand(N)*2 - 1.0

ratio = np.mean((x**2 + y**2 < 1.0))

print(N, ratio*4)

원의 면적을 구해보는 위의 예처럼, 일반적으로 몬테카를로 방법은 마구잡이 변수를 이용한 수치 적분이라고 생각할 수 있다. 예를 들어, [eqn: int_0^5 e ^{-x ^{2}} dx]의 값을 몬테카를로 방법을 이용해 구해보는 문제를 생각해보자. 위의 코드와 같은 방법으로 [eqn: x in [0,5] , y in [0,1]]의 직사각형 영역에서 마구잡이로 [eqn: (x,y)]를 여러 번 택해 이중 [eqn: y < e^-x^2]을 만족하는 점의 숫자를 세어서 적분값을 구할 수도 있다. 하지만, 조금만 생각해보면 이 방법은 심각한 결함이 있다는 것을 알 수 있다. [eqn: e^-x^2]이 급격히 줄어드는 꼴이라서 [eqn: x]가 큰 영역에서는 이 부등식을 만족하는 [eqn: (x,y)]가 거의 택해지지 않기 때문이다.

== 예제 1 ==

위에서 설명한 방법을 따라 [eqn: (x,y) in [0,5] times[0,1]]를 만족하는 무작위 수 [eqn: (x,y)]를 [eqn: N]번 택해서 [eqn: int_0^5 e ^{-x ^{2}} dx]의 값을 계산하는 코드를 작성하시오.

좀 더 나은 다른 방법도 있다. [eqn: x in[0,5]]를 만족하는 [eqn: x]를 [eqn: N]개 무작위로 택하면, 인접한 두 [eqn: x]사이의 간격은 [eqn: Delta x approx 5/N]이 된다. 이를 이용하면 적분 [eqn: int_0^5 e ^{-x ^{2}} dx]의 근사값은 [eqn: 5 over N sum_{i=1}^N e^-x_i^2]로 어림할 수 있다. 이 방법으로 적분값을 구하는 프로그램은 아래와 같다. 프로그램을 실행해서 위 예제의 결과와 비교해보라.

import numpy as np, matplotlib.pyplot as plt

for N in 10**np.array([1,2,3,4,5,6,7]):

x = np.random.rand(N)*5

y = np.exp(-x**2)

print(N, np.sum(y)*5.0/N)

중요성 샘플링은 중요한 부분에서 더 촘촘하게 샘플링을 하는 방법이다. 중요성 샘플링(오른쪽)을 하면 균일한 샘플링(왼쪽)보다 더 정확한 적분값을 얻게 된다.

다음에는 몬테카를로 방법에서 널리 쓰이는 중요성 샘플링(importance sampling)에 대해 알아보자. 위에서 생각한 함수 [eqn: e^-x^2]은 [eqn: x]가 커지면서 급격히 그 값이 줄어들게 되므로, 최종 적분값에 대한 대부분의 기여는 [eqn: x]가 작을 때에 주로 생기게 된다. 이처럼 함수가 특정 영역에서만 큰 값을 갖는 경우에는 이 영역에서는 좀 더 촘촘하게 샘플링을 하고, 영역 밖에서는 성기게 샘플링을 하는 것이 보다 정확한 수치 적분을 얻기에 유리하다는 것을 알 수 있다. 이처럼, 영역의 중요성에 따라 샘플링을 달리하자는 것이 바로 중요성 샘플링의 기본 아이디어이다. 계산하고자 하는 함수 [eqn: f(x)]의 구간 [eqn: [a,b]]에서의 적분 [eqn: I = int_a^b f(x) dx]을 위의 간단한 프로그램에서는 아래와 같이 어림해서 계산했다.

[eqn: I APPROX {b-a} over {N} sum _{i=1} ^{N} f(x _{i} )]

만약 확률분포함수 [eqn: p(x)]를 이용한다면 [eqn: x]를 샘플링하면, 구간 [eqn: [x, x+Delta x ]]에서 [eqn: x]가 택해질 확률은 [eqn: p(x) DELTA x]이고, 모두 [eqn: N]개의 [eqn: x]를 샘플링하므로 구간 [eqn: [x, x+Delta x ]]안에서 샘플링한 [eqn: x]의 개수는 [eqn: n = Np(x) Deltax]가 된다. 이 구간 안에서 샘플링된 [eqn: x_1, x_2 , cdots , x_n]을 생각하면, 구하고자 하는 적분값 중 구간 [eqn: [x, x+Delta x]]이 기여하는 부분은 [eqn: f(x) DELTA x APPROX left[ 1 overn sum_{i=1}^n f(x_i ) right] Deltax = 1 over{Np(x)} sum_{i=1}^n f(x_i ) approx 1 over {N} sum_{i=1}^n {f(x_i )}over{p(x_i )]로 적을 수 있다. 이를 전체 구간 [eqn: [a,b]]에 대해 모두 더하면 아래의 식을 얻게 된다.

[eqn: I approx 1 over N sum_{i=1}^N {f(x_i )} over {p(x_i ) }]

다음에는 [eqn: I = int _{0} ^{5} {} e ^{-x ^{2}} dx]의 값을 중요성 샘플링을 이용해 계산해보자. 중요성 샘플링을 이용하기에 가장 좋은 함수는 물론 [eqn: p(x) = A e ^{-x ^{2}}]임은 분명하지만, 비례상수 [eqn: A]의 값이 다름 아닌 [eqn: I]이므로 [eqn: p(x) = A e ^{-x ^{2}}]를 이용할 수는 없다. 아래에서는 대신 [eqn: p(x) = Ae^-x]를 이용하고자 한다. 상수 [eqn: A]는 [eqn: int_0^5 p(x) dx = A (1-e^-5 )= 1]이므로 [eqn: A = 1 over{1-e^-5 }]이다. 위에서 얻은 식 [eqn: I approx 1 over N sum_{i=1}^N {f(x_i )} over {p(x_i ) }]을 이용하기 위해서는 먼저 [eqn: x]를  확률분포함수 [eqn: p(x) = Ae^-x]를 이용해 샘플링 해야 한다.

임의의 확률분포 [eqn: p(x)]를 갖는 마구잡이 수를 얻는 간단한 방법이 있다. [eqn: y]가 구간 [0,1]에서 일정한 값을 갖는 균일한 확률분포 [eqn: p(y)=1]를 갖는 확률변수라면 [eqn: p(x)dx=p(y)dy`=`dy]이므로 [eqn: y`= int _{0} ^{x} {} p(x)dx= {A} int _{0} ^{x} {} e ^{-x} dx= {A} (1-e ^{-x} )]를 얻고, 이를 정리하면 [eqn: 1- y over A =e ^{-x}]이므로 [eqn: x=-ln LEFT ( 1`- y over A right )]이다. 즉, 0과 1사이의 균일한 분포를 가진 마구잡이 수 [eqn: y]를 생성하고 이에 해당하는 값 [eqn: x = -ln left( 1 - y over A right)]를 계산하고는 이를 이용해 [eqn: I APPROX {1} over {N} sum _{i=1} ^{N} {f(x _{i} )} over {p(x _{i} )} = 1 over N sum_{i=1}^N {e^{-x_i^2 } } over {A e^{-x_i }]를 계산하면 [eqn: I = int _{0} ^{5} {} e ^{-x ^{2}} dx]의 어림값을 얻게 된다. 아래의 프로그램을 살펴보라.

import numpy as np, matplotlib.pyplot as plt

A = 1.0/(1.0 - np.exp(-5.0))

for N in 10**np.array([1,2,3,4,5,6,7]):

y = np.random.rand(N)

x = -np.log( 1.0 - y/A)

I = np.sum( np.exp(-x*x)/np.exp(-x) )/A/N

print(N, I)

프로그램이 출력한 값을 [eqn: I= int _{0} ^{5} {} e ^{-x ^{2}} dx = 1 over 2 sqrt{pi} rmerf(5) approx 0.886227]과 비교해보라.

== 예제 2 ==

[eqn: I = int _{0} ^{5} {} e ^{-x ^{2}} dx]의 값을 확률분포 [eqn: p(x)`=`1`-`x/5]를 이용한 중요성 샘플링을 통해 구하시오.

통계역학에서의 몬테카를로 방법의 이용

통계역학의 바른틀 모듬(canonical ensemble)에서 계가 특정 상태 [eqn: vec x]에 있을 확률은 이 상태에 해당하는 해밀토니안 [eqn: H(vec x )]를 이용해 볼츠만 확률분포 [eqn: {e ^{- beta H( {vec{x}} )}} over {Z}]로 주어진다. 여기서 [eqn: Z = int d vec x e^{-beta H(vec x ) }]는 분배함수라 불린다. 통계역학을 이용해 물리량 [eqn: O(vec x )]의 바른틀 모듬의 평균을 구하는 식은  다음과 같이 적분의 형태로 적힌다.

[eqn: <O( {vec{x}} )>`=` {1} over {Z} int _{} ^{} {} d {vec{x}} O( {vec{x}} )e ^{- beta H( {vec{x}} )}]

즉, 물리량의 열역학적인 평균값을 구하는 것은 수학적으로는 적분을 수행하는 과정일 뿐이다. 위의 식에서 적분이 수행되는 벡터공간은 거시적인 물리계의 경우에 엄청나게 높은 차원을 가지게 되어서, 단순한 방법으로는 적분을 하는 것이 불가능하다. 예를 들어, 3차원에서 움직이는 100개의 입자가 있는 경우, 입자하나당 위치와 운동량이 각각 3차원 벡터이므로, 위의 적분은 모두 600차원 벡터공간에서의 적분이 된다. [eqn: vec x]가 연속적인 값이 아닌, 전자의 스핀처럼  +1, -1의 두값만을  가지는 경우에는 적분이 아니라 합계를 구하는 방법을 이용하게 되어,   [eqn: <O( vec sigma )> = 1 over Z sum_{sigma_1 = plusminus 1 } sum_{sigma_1 = plusminus 1 } cdots sum_{sigma_100 = plusminus 1 } O( vec sigma ) e^{-H( vec sigma ) }]를 이용할 수 있다. 하지만, 이 경우에도 더해야 하는 항의 개수가 [eqn: 2^100]이어서 컴퓨터를 이용하더라도 정확한 값을 계산하는 것은 불가능에 가깝다.

중요성 샘플링을 이용한 몬테카를로 방법을 적용하면 위의 난점을 해결할 수 있다. 거시적인 계의 경우에는 [eqn: vec x]가 존재하는 높은 차원의 벡터공간의 영역 안에서 [eqn: < O( vecx ) >]에 유의미하게 기여하는 영역은 극히 작기 때문이다. 이 부분에서 주로 샘플링이 이루어지도록 중요성 샘플링을 적용하면 아주 효율적으로 열역학적인 양들을 계산할 수 있게 된다.

적분하고자 하는 함수 [eqn: f(x)]에 비례하도록 중요성 샘플링의 확률분포를 택하는 것이 유리하다고 할 수 있지만, 실제로 이를 이용할 수는 없다. 즉, [eqn: p(x) = A f(x)]로 적으면, 비례상수 [eqn: A]는 [eqn: int _{} ^{} {} p(x)dx=A int _{} ^{} {} f(x)dx=1]의 조건을 이용해 얻어야 하는데, [eqn: 1/A`= int f(x) dx]이므로, 구하고자 하는 적분값 [eqn: int f(x) dx]를 알아야만 [eqn: p(x) = A f(x)]를 이용할 수 있기 때문이다. 통계역학의 몬테카를로 계산에서 가장 널리 쓰이는 방법인 메트로폴리스(Metropolis) 알고리즘은 바로 이 측면에서 신기한 방법이다. [eqn: A]의 값을 구하지 않고도 [eqn: p(x) = A f(x)]를 이용할 수 있다. 바른틀 모듬에서의 열역학적인 평균에 대해 바꿔 얘기하면,  [eqn: <O( {vec{x}} ) > = {1} over {Z} int _{} ^{} {} d {vec{x}} O( {vec{x}} )e ^{- beta H( {vec{x}} )}]에서 [eqn: Z]를 모르고도 [eqn: p( vec x ) = e^{-beta H( vec x ) } over Z]를 이용해 중요성 샘플링을 하는 놀라운 방법이다. 중요성 샘플링의 확률분포를 이처럼 볼츠만 확률분포로 택하면 평균값은 아래의 식으로 아주 간단한 산술평균의 꼴로 적히게 된다.

[eqn: <O( {vec{x}} )> = {1} over {Z} int _{} ^{} {} d {vec{x}} O( {vec{x}} )e ^{- beta H( {vec{x}} )} approx 1 over ZN sum_i=1^N O(vec x_i ) e^{-beta H(vec x _i ) } over {p ( vec x _i ) } = 1 over N sum_{i=1}^N O( vec x _i )]

메트로폴리스 알고리즘

시스템이 [eqn: mu]라는 상태에 있을 확률을 [eqn: P_mu]라고 적자. 그리고 시스템이 [eqn: mu]에서 다른 상태 [eqn: nu]로 전이할 확률(전이확률, transition probability)을 [eqn: W( mu ` -> ` nu )]라고 하자. [eqn: P_mu]가 시간([eqn: t])에 따라 변해가는 것을 기술하는 것이 바로 아래의 으뜸방정식(master equation)이다.

[eqn: {dP _{mu }} over {dt} = sum _{nu} ^{} P _{nu } W(` nu ` -> ` mu )- sum _{nu } ^{} P _{mu } W( mu ` -> ` nu )]

상태 [eqn: mu]에 있을 확률은 다른 상태 [eqn: nu]가 [eqn: mu]로 전이하는 과정에서는 커지고, [eqn: mu]가 [eqn: nu]로 전이하는 과정에서는 줄어드는 것을 생각하면 으뜸방정식의 의미를 쉽게 이해할 수 있다. 메트로폴리스 알고리즘은 시간이 흐르면서 확률분포 [eqn: P_mu]가 평형상태의 확률분포인 [eqn: P _{mu } ^{eq} = e^{-beta H_mu } over Z]를 향해 접근하도록 하는 방법이다. 평형상태에 도달한 다음에는 [eqn: P_mu]는 시간에 따라 더 이상 변하지 않으므로 위의 으뜸방정식에서 [eqn: dP_mu^eq over dt = 0]이 되고, 따라서 아래의 식을 만족하게 된다.

[eqn: sum _{nu} ^{} LEFT [ P _{nu } ^{eq} W( nu ` -> ` mu )-P _{mu } ^{eq} W( mu ` -> ` nu ) RIGHT ] =0].

다음에는 위의 조건이 전체의 총합이 아닌 개개의 항 하나하나에도 성립하다는 조건(detailed balance라고 부른다)을 추가로 요구해 다음의 조건을 얻는다.

[eqn: {W( nu rightarrow mu )} over {W( mu rightarrow nu )} = {P_mu^eq} over {P_nu^eq } = e^{ -beta (H_mu - H_nu ) }]

즉, [eqn: {W( nu -> mu )} over {W( mu -> nu )} =e ^{- beta (H _{mu } -H _{nu } )}]를 만족하도록 전이확률 [eqn: W(mu rightarrow nu )]을 택하면, 위의 으뜸방정식을 통해 확률분포 [eqn: P_mu]는 시간이 지나면서 [eqn: P_mu^eq = e^{-beta H_mu } over Z]로 저절로 다가서게 되고, 이후에는 중요성 샘플링이 바로 이 볼츠만 확률분포를 따라서 이루어지게 된다. 시스템의 현재 상태를 [eqn: mu]라고 하고 다음 전이할 상태를 [eqn: nu]라고 하면 두 상태의 에너지 차이는 [eqn: Delta E = H_nu - H_mu]가 되는데, 표준적인 메트로폴리스 알고리즘에서는 아래의 전이확률을 이용한다(위의 detailed balance 조건을 만족함을 확인해보라).

만약 [eqn: Delta E leq 0]이면, [eqn: W( mu rightarrow nu ) = 1]

만약 [eqn: Delta E > 0]이면, [eqn: W(mu rightarrow nu ) = e^{-beta Delta E].

메트로폴리스 알고리즘을 이용한 몬테카를로 방법의 적용: 이징 모형

이징 모형은 통계물리학의 발전에서 큰 역할을 한 간단한 모형이다. 이징 모형의 해밀토니안은 아래의 형태로 적힌다.

[eqn: H= - J over 2 sum _{ij} ^{} A _{ij`} S _{i} S _{j}]

[eqn: i]번째 스핀값 [eqn: S_i = plusminus1]의 두 값을 가지며, 강자성 상호작용을 하는 경우 [eqn: J]는 0보다 크다. 상호작용의 구조를 기술하는 행렬(인접 행렬, adjacency matrix)의 성분 [eqn: A_ij]는 두 스핀 [eqn: S_i]와 [eqn: S_j]가 상호작용하고 있는 경우에는 1의 값을 그렇지 않은 경우에는 0의 값을 가진다. 두 스핀의 상호작용은 반대 방향으로도 있으므로 [eqn: A_ij]는 대칭행렬이고 두 스핀사이의 상호작용이 합계를 계산할 때 두 번씩 등장하므로 전체를 2로 나눴다. 2차원의 사각격자의 형태로 놓여 가장 가까운 스핀하고만 상호작용을 하는 최인접 상호작용 이징 모형 뿐 아니라 임의의 상호작용의 구조를 인접행렬 [eqn: A_ij]를 이용해 구현할 수 있다. 이징 모형은 특정한 온도에서 상전이를 보여주는데, 온도가 상전이온도보다 낮은 경우에는 거시적인 숫자의 스핀들이 같은 방향을 가리켜서 [eqn: m = 1 overN sum_i S_i]의 값이 0보다 큰 유한한 값을 갖고, 높은 온도에서는 스핀들의 방향이 뒤죽박죽이 되어서 [eqn: m]의 값이 0이 된다. 이처럼 이징 모형의 상전이는 질서도라 불리는 [eqn: m]의 평균값을 이용해 연구할 수 있다. 몬테카를로 방법으로 이징 모형의 상전이를 알아보는 프로그램을 작성해보자. 메트로폴리스 알고리즘에서는 현재의 주어진 스핀 상태 [eqn: {vec{S}} = LEFT { S _{1} ,S _{2} ,S _{3} , CDOTS S _{N} RIGHT }]로부터 다른 상태 [eqn: {vec{S}} prime = LEFT { S _{1}prime ,S _{2}prime ,S _{3}prime , CDOTS S _{N}prime RIGHT }]로의 전이확률로 만약 [eqn: DELTA E = H( {vec{S}} ‘)-H( {vec{S}} )]의 값이 0보다 작으면 [eqn: W( vec S rightarrow vec S ‘ ) = 1], 0보다 크면 [eqn: W( {vec{S}} -> {vec{S}} ‘)= e^{-beta Delta E }]를 이용한다. [eqn: vec S]로부터 [eqn: vec S prime]을 만드는 방법으로는 다양한 방법이 있지만, 가장 단순한 경우로서 [eqn: N]개의 스핀 중 딱 하나의 스핀만이 뒤집힌 상황을 생각하자. 즉, [eqn: vec S prime]에 들어있는 스핀 중 [eqn: k]번째 스핀 하나만 [eqn: S _{k} prime =-S _{k}]를 만족한다.

이징 모형의 몬테카를로 계산은 아래의 절차를 따라 수행된다.

전체 [eqn: N]개의 스핀 중 임의로 하나를 골라([eqn: k]번째 스핀이라 하자), [eqn: S_k]를 [eqn: -S_k]로 바꾸자. 이렇게 만들어진 스핀 상태를 [eqn: vec S prime]이라 하자.

에너지의 차이 [eqn: Delta E = H(vec S prime ) - H( vec S )]를 계산한다.

만약 [eqn: Delta E leq 0]이면 [eqn: vec S]에서 [eqn: k]번째 스핀을 뒤집는다. 즉, [eqn: S _{k} ` -> `-S _{k}].
만약 [eqn: Delta E > 0]이면, [eqn: e^{-beta Delta E}]의 확률로 [eqn: vec S]에서 [eqn: k]번째 스핀을 뒤집는다. 즉, [eqn: S _{k} ` -> `-S _{k}].
(0과 1사이의 마구잡이 수 [eqn: r]을 생성하고는 [eqn: r < e^{-beta Delta E]을 만족하면 스핀을 뒤집으면 된다.)

질서도 [eqn: m = 1 over N sum_i S_i]를 계산해 저장한다.

다시 단계 1로 이동하고 전체 과정을 여러 번 반복한다.

모든 단계가 끝난 후 질서도의 평균값 [eqn: < m >]을 출력한다.

위의 전 과정을 온도 [eqn: T]를 바꿔가면서 반복한다.

전역적인 상호작용을 하는 이징 모형은 모든 [eqn: i, j]에 대하여 [eqn: A_ij = 1 over N]를 이용해 구현할 수 있다. 1이 아니라 [eqn: 1/N]을 이용한 이유는, 시스템의 해밀토니안이 [eqn: N]에 비례해야하기 때문이다. 아래의 프로그램은 온도 [eqn: T]의 함수로 [eqn: m]과 비열을 구해 그래프로 그린다.

import numpy as np, matplotlib.pyplot as plt

N = 1000 # number of spins

MCstep = 10000

def calE(A,S):

return(-0.5*S@A@S.T)

def MCrun(A, S, N, T, E):

Snew = np.copy(S)

k = np.random.randint(N) # Pick a spin to update

Snew[k] *= -1

Enew = calE(A,Snew)

DeltaE = Enew - E

if (DeltaE <= 0.0): S[k] *= -1; E = Enew

elif (np.random.rand() < np.exp(-DeltaE/T)): S[k] *= -1; E = Enew

return(E)

A = np.ones([N,N])/N # adjacency matrix for global coupling

S = 2*(np.random.rand(N) > 0.5) -1 # initialization of spin vector

xarr = []; yarr1 = []; yarr2 = []

E = calE(A,S)

for T in np.arange(0.1,2.0,0.1): # temperature loop begins

# Equilibration

for n in range(MCstep):  # MC loop begins

E = MCrun(A,S,N,T,E)

# Measurement

m = 0.0; sumE = 0.0; sumE2 = 0.0

for n in range(MCstep):  # MC loop begins

E = MCrun(A,S,N,T,E)

m += np.abs(np.mean(S))

sumE += E; sumE2 += E*E

print(T, m/MCstep)

sumE /= MCstep

sumE2 /= MCstep

xarr.append(T)

yarr1.append(m/MCstep)

yarr2.append((sumE2 - sumE*sumE)/T/T/N)

plt.plot(xarr,yarr1)

plt.show()

plt.plot(xarr,yarr2)

plt.show()

== 과제 1 ==

전역적인 상호작용을 하는 이징 모형의 경우 [eqn: M = sum_i S_i]를 변수로 이용하면 [eqn: H = - J over 2N M^2]이 된다. 이식을 이용해 위의 프로그램을 보다 효율적으로 만들어보시오. [eqn: k]번째 스핀을 뒤집는 시도와 관련된 에너지차이 [eqn: DeltaE]도 [eqn: M]을 이용하면 쉽게 계산할 수 있음도 이용하시오.

== 과제 2 ==

이차원 사각격자 위에서 최인접 스핀하고만 상호작용을 하는 이징 모형을 생각하자. 주기적인 경계조건이 적용된 [eqn: A_ij]를 구현하고, 몬테카를로 계산을 통해 [eqn: m]을 얻어 온도의 함수로 그래프로 그리시오. 상전이 온도를 추정하시오.